---

pca_df_total = read.csv('data/total_pre_pca.csv', header = TRUE) #climate values from all 19 BIOCLIM variables at both Xant and Yucc unique locations

pca_clim_total = prcomp(pca_df_total, center = T, scale. = T) # PCA dimension analysis for 19 bioclim variables and the elevation from all locations of the night lizard and the joshua tree 

# visualize the results

# screeplot
fviz_eig(pca_clim_total) # we see that the first two component contain a substantial amount of the total variance

# visualizing variables
fviz_pca_var(pca_clim_total,
             col.var = "contrib", 
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = T) #We see that precip and temp variables do split between pca directions, but not for all cases

# Here, we upload the appended dataframe where both night lizard and joshua tree features were rotated into the above PCA
total_pca_complete = read.csv('data/total_pca.csv', header = TRUE) # This is the final dataframe of PCA components

# some of this code, done in a seperate markdown file: 

# Once we had the PCA from the total habitat, we rotated each species's set of BIOCLIM features with the total PCA from above. These steps are ommitted here for brevity, and a susbet of the code for this is commented out below:
# yucca = scale(yuc_climate, pca_clim_total$center, pca_clim_total$scale) %*% pca_clim_total$rotation
# xan = scale(xan_climate, pca_clim_total$center, pca_clim_total$scale) %*% pca_clim_total$rotation
# write.csv(xan, file = 'xanpcatotal.csv') #write.csv(yucca, file = 'yucpcatotal.csv')

Visualizing results

p_total <- plot_ly(total_pca_complete, x = ~PC1, y = ~PC2, z = ~PC3, color = ~X, colors = c('#BF382A', '#0C4B8E')) %>%
  add_markers() %>%
  layout(scene = list(xaxis = list(title = 'PC1'),
                     yaxis = list(title = 'PC2'),
                     zaxis = list(title = 'PC3')))
p_total

xan_set1 <-
  total_pca_complete %>%
  filter(X == 0) 
xan_set1 = xan_set1[,2:4]
yucca_set1 <-
  total_pca_complete %>%
  filter(X == 1) 

yucca_set1 = yucca_set1[,2:4]

hull_xan1 = convhulln(xan_set1, output.options = 'FA')
hull_xan1$vol #vol = 531.6

## [1] 531.5689

hull_yuc1 = convhulln(yucca_set1, output.options = 'FA')
hull_yuc1$vol #vol = 516.9

## [1] 516.8733

xanmesh = to.mesh3d(hull_xan1)
yucmesh = to.mesh3d(hull_yuc1)

wire3d(xanmesh, meshColor = c("faces"), color = "red")
dot3d(xanmesh, meshColor = c("faces"), color = "red")
wire3d(yucmesh, meshColor = c("faces"), color = "blue")
dot3d(yucmesh, meshColor = c("faces"), color = "blue")